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Abstract. We present a numerical method for the approximation of solutions for the class of stochastic differential equations driven by 
Brownian motions which induce stochastic variation in fixed directions. This class of equations arises naturally in the study of population 
processes and chemical reaction kinetics. We show that the method constructs paths that are second order accurate in the weak sense. The 
method is simpler than many second order methods in that it neither requires the construction of iterated Ito integrals nor the evaluation of 
any derivatives. The method consists of two steps. In the first an explicit Euler step is used to take a fractional step. The resulting fractional 
I point is then combined with the initial point to obtain a higher order, trapezoidal like, approximation. The higher order of accuracy stems 

from the fact that both the drift and the quadratic variation of the underlying SDE aie approximated to second order. 

(N 

^ 1. Introduction We consider the problem of constructing accurate approximations on bounded time 

^ intervals to solutions of the following family of stochastic differential equations (SDEs) 

^ M 

^ dX{t)^b{X{t))dt + Y,'^k{X{t))i^kdWkit), 

k=l ^ 

where 6: M"* — > M'^, ak-R'^ — > R>q, Vk G IR'*, and Wk{t) are one-dimensional Wiener processes. Thus, 
randomness is entering the system in fixed directions Vk, but at variable rates ak{X{t)). Precise regularity 
"j^ conditions on the coefficients will be presented with our main results in Section|2] 

^ The algorithm developed in this paper is a trapezoidal-type method and consists of two steps; in the first 

an explicit Euler step is used to take a fractional step and in the second the resulting fractional point is used 
in combination with the initial point to obtain a higher order, trapezoidal like, approximation. We will prove 
that the method developed is second order accurate in the weak sense. Because the method developed here 
produces single paths, it is natural to allow variable step-sizes; this is in contrast to Richardson extrapolation 
f — techniques ( ll20l ). Finally, it is important to note that while the method presented in this paper is applicable to 

only a sub-class of SDEs, that sub-class does include systems whose diffusion terms do not commute, which is 
a classical simplifying assumption to obtain higher order methods (See ll9l [T4l '). 
\Q The method we propose is in some sense similar to the classical predictor-corrector There have already 

been a number of such methods proposed in the stochastic context to produce higher-order methods (see fTsl 
2^ [mia). In a general way, all of these methods require the simulation of iterated Ito integrals and sometimes 

IT. need derivatives of the diffusion terms. If one only cares about weak accuracy, it is possible to use random 

^ variables which make these calculations easier and computationally cheaper That being said, the complexity 

and cost of such calculations is one of the main impediments to their wider use. By assuming a certain structure 
for we are able to develop a numerical method which we hope is more easily applied and implemented. 

Though a specific structure of is assumed, it is a structure which arises naturally in a number of 
settings. For example, our method will be applicable whenever d — 1. Also, we note that diffusion approxi- 
mations to continuous time Markov chain models of population processes, including (bio)chemical processes, 
satisfy \\.\\ . As stochastic models of biochemical reaction systems, and, in particular, gene regulatory systems, 
are becoming more prevalent in the science literature, developing algorithms that utilize the specific structure of 
such models has increased importance (UlEl). Furthermore, in Section[8] we quote a result from the literature 
which states that any system with uniformly elliptic diffusion can be put in the form of without changing 
its distribution. 
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The topic of this paper is a method that produces a weak approximation rather than a strong approximation 
in that the approximate trajectory is produced without reference to an underiining Wiener process trajectory. 
We see this as an advantage. Except for appHcations such as fihering or certain problems of collective motion 
for stochastic flows, one is usually simply interested in generating an accurate draw from the distribution on 
C{[0,T],R'^) induced by ( fTT) . This is different than accurately reproducing the Ito map W H> X{t,W) 
implied by The second is referred to as strong approximation. In our opinion such approximations are 

usually unnecessary and lead to a concept of accuracy which is unnecessarily restrictive. In |7|, it is discussed 
that without accurately estimating second order Ito integrals one cannot produce a strong method of order 
greater then 1/2. If the vector fields commute, then this restriction does not apply and higher order strong 
methods are possible. While the term "strong approximation" is quite specific, the term "weak approximation" 
is used for a number of concepts. Here we mean that the joint distribution of the numerical method at a fixed 
number of time points converges to the true marginal distribution as the numerical grid converges to zero. If 
this error goes to zero as the numerical mesh size to the power p in some norm on measure then we say the 
method is of order p. This should be contrasted with talking about the rate at which a given function of the path 
converges. 

The outline of the paper is as follows. In Section|2]we present our algorithm together with our main results 
concerning its weak error properties. In Section[3]we give the intuition as to why the method should work. In 
Section|4]we give the delayed proof of the local error estimates for the method which were stated in Section|2] 
In Section |5] we provide examples illustrating the performance of the proposed algorithm. In Section |6] we 
discuss the effect of varying the size of the first fractional step of the algorithm. In SectionlTlwe compare one 
step of the algorithm to one step in a Richardson extrapolation type algorithm. In Section |8] we show how, at 
least theoretically, the method can be applied to any uniformly elliptic SDE. Finally, an appendix contains a 
tedious calculation needed in Section!?] 

2. The numerical method and main results Throughout the paper, we let X{t) denote the solution to 
( |l.l| l and Yi denote the computed approximation at the time ti for the time discretization = to < < • • • • We 
begin both from the same initial condition, namely X{0) = lo = xq. Let {jyj^; ' ^ ^ {1, . . . , A/}, « G N} 
be a collection of mutually independent Gaussian random variables with mean zero and variance one. It is 
notationally convenient to define [x]^ = xV = max{a;, 0}. 

We propose the following algorithm to approximate the solutions of 

Algorithm. ( Weak 0-Midpoint Trapezoidal ) Fixing a 6* e (0, 1), we define 

Next fixing a discretization step h, for each i G {1, 2, 3, . . . } we repeat the following steps in which we first 
compute a 9-midpoint y* and then the new value Yi. 

M 

Step 1. y* = + biY,^i)9h + ^ ak{Y,^i) Vk v[2^ 

k = l 

M 



Step 2. Y, = y* + {aib{y*) - a2b{Y,_i)){l - d)h + J] V ["i'^fc(2^*) - {I - 6)h. 

fc=i 

Remark 2.1. Notice that on the ith-step y* is the standard Euler approximation to X{9h + (i — l)h) starting 
from Yi^i at time (i — l)h M3\l . 

Remark 2.2. Notice that for all 9 € (0, 1) one has ai > a2 and a\ — a2 = 1- It is reasonable to ask which 
9 is best. Notice that when 9 — \/2 both ai and ai are minimized with values ai — 2 and ai — 1. This likely 
has positive stability implications. From the point of view of accuracy 9 = 1/2 also seems like a reasonable 
choice as it provides a central point for building a balanced trapezoidal approximation, as will be explained in 
Section^ Further, picking a 9 close to 1 or increases the likelihood that the term [aiaf{y*) — Q!20'fc(^i-i)]^ 

2 



will be zero, which will lower the accuracy of the method. If instability due to stiffness is a concern, one might 
consider a 6 closer to one as that would likely give better stability properties being closer to an implicit method. 
In general, 6=1/2 seems like a reasonable compromise, though this question requires further investigation 
and will be briefly revisited in Section [6] 

For simplicity, we will restrict ourselves to the case when b and the ak are in C^(M''), the space of bounded 
functions whose first through sixth derivatives are continuous and bounded. In general, we will denote by 
C'^(M'^) the space of bounded, continuous functions whose first k derivatives are bounded and continuous. For 
/ e C'' {R'^), we define the standard norm 

ll/lk = sup{|/(x)|, \da,fix)\ ■.xeR'^,a = {ai, . . . , aj), a, G {!,..., < k} . 

It is notationally convenient to define the Markov semigroup Vt'- associated with \\.\\ by 

{Vtf){xY^^J{X{t)) (2.2) 

where ^(0) = x and Markov semigroup Ph'. — > associated with a single full step of size h of the 
numerical method by 

where Iq = V- Clearly ||P/i/||o < ll/llo ™d ||Ph/||o < ll/llo- It is also a standard fact, which we summarize 
in Appendix [b] that in our setting for any t > and k e N if 6, cti , . . . gm G then there exists a C = 
C(r, fc, 6, a) so that ||7^f/||A.. < C||/||fc is true for all t < T. All of these can be rewritten succinctly in the 
induced operator norm from C'' — > C'' as ||7't||fc->.fc < C, ||7't||o->.o < 1 and ||P/i||o^o < 1- Analogously, for 
any linear operator L: C'' ^ we will denote the induced operator norm from C'' — > by which 
is defined by 

= sup 



/sc'^Jt^o ll/lk- 

The following two theorems are the principle results of this article. They give respectively the weak local 
and global error of the Weak Trapezoidal method. 

Theorem 2.3 (One-step approximation). Assume that b G and for all k, ct^ e with vai^ crk[x) > 
0. Then there exists a constant K so that 

< Kh^ (2.3) 

for all h sufficiently small. 

From this one-step error bound, it is relatively straight-forward to obtain a global error bound. The follow- 
ing result shows that our approximation scheme gives a weak approximation of second order. 

Theorem 2.4 (Global approximation). Assume that b £ and for all k, ak G with inix <Tk{x) > 0. 
Then for any T > there exists a constant C{T) such that 

sup WVnh - Phk^o < C{T)h^. (2.4) 

0<n<T/h 



Proof. We begin by observing that 

n 
fc=l 
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and hence since supQ^^^j. ||7's||6-i.6 < C{T) and ||P;J'||o->.o < 1> using \23) we have that for any n with 

< n < T/h 

n 

\\Vnh - Phk^O < \\Pt^\\0^0\\Ph - Vh\\6^0\\Vh(n-k)\\6^6 

fc=l 
n 

< C{T)Kh^ = KTC{T)h^ = C(r)/i2. 

k=l 

□ 

Remark 2.5. r/ie restriction that 'mix c/c(x) > can likely be relaxed if one has some control of the behavior 
of the solution around the degeneracies of ak{x). This assumption is made to keep the proof simple with easily 
stated assumptions. 

3. Why the method works We now give two different, but related, explanations as to why the Weak 
6'-Midpoint Trapezoidal Algorithm is second order accurate in the weak sense. 

3.1. A first point of view Inserting the expression for y* from Step 1 of the Weak ^-Midpoint Trape- 
zoidal Algorithm into Step 2 and disregarding the diffusion terms yields 



(3.1) 



Considering ( |3.1| l, we see that when w 1 we recover the standard theta method (not to be confused with our 
use of 9) with theta =1/2, which is known to be a second order method for deterministic systems. When 6 = 
1/2, we recover the standard trapezoidal or midpoint method. For 9 ^ 1/2, we simply have a trapezoidal rule 
where a fractional point of the interval is used in the construction of the trapezoid. We will argue heuristically 
that the Weak Trapezoidal Algorithm handles the diffusion terms similarly. We also note that p.2[ ) shows that 
our algorithm can be understood as an approximation to the two-step Taylor series where 6* is a parameter used 



to approximate the second derivative. This idea will be revisited in the proof of Theorem 2.3 
Equation ( |l.l| i is distributionally equivalent to 

X{t)=X{Q)+ b{X{s))ds + Yuk \^,,2,x{s))){u)Yk{duy.d.s), (3.3) 

Jo j,^^ Jo JO 

where the Yfe are independent space-time white noise processed and all other notation is as before, in that 
solutions to p.3| l are Markov processes that solve the same martingale problem as solutions to \\.\\ \ that is, 
they have the same generator (|6|). In order to approximate the diffusion term in p.3| l over the interval [0, h), 
we must approximate yA;(^[o,h) i^t)) where ^[o,/i) (c^) is the region under the curve a1.{X{t)) for < t < h. 

We consider a natural way to approximate X{h) and focus on the double integral in ( |3.3[ ) for a single 
k. We also take 9 — 1/2 for simplicity and simply note that the case 9^1/2 follows similarly. We begin 
by approximating the value X{h/2) by y* obtained via an Euler approximation of the system on the interval 
[0, h/2). To do so, we hold X{t) fixed at X{Q) and see that we need to calculate Yfe (Region 1), where Region 
1 is the grey shaded region in Figure [3T{ a). Because 

yfe(Regionl) ^ iV(0,afe2(X(0))V2) = afe(X(0))./| iV(0, 1), 



'More precisely, the are random measures on [0, oo)^ such that if A,B (Z [0, oo)^ with Ar\ B = % then Yk{A) and Yj.(B) are 
each independent, mean zero Gaussian random variables with variances Area(A) and Area(i3), respectively. Integration with respect to 
this field can be defined in the standard way beginning with adapted simple functions which are fixed random variables on fixed rectangular 
sets and then extending by linearity after the appropriate Ito isometry is established. 
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Region 5 



^ h_ h_ M h_ — ►.I h_ — ^ h_ — h_ — 

2 2 2 2 2 2 

(a) First step (b) Desired second step (c) Used second step 

Fig. 3.1: A graphical depiction of the Weak Trapezoidal Algorithm with 9 = 1/2. In (a) the region of space- 
time used in the first step of the Weak Trapezoidal Algorithm is depicted by the grey shaded Region 1. In 
(6) the desired region to use, in order to perform a trapezoidal approximation, would be Region 2. However 
we have used Region 3 in our previous calculation and this is analytically problematic to undo. In (c), where 
V = crliu*) — a1{X{0)), we see that Region 5 gives the correct amount of new area wanted as subtracting off 
the area of Region 4 "offsets" the used area of Region 3. The case 9 ^ 1/2 is similar. 



we see that this step is equivalent in distribution to Step 1 of Algorithm |l.l| (Here and in the sequel, " = " 
denotes "equal in distribution.") 

If we were trying to determine the area under the curve (jl{X{t)) using an estimated midpoint y* for a 
deterministic X{t), one natural (and common) way would be to use the area of Region 2, where Region 2 is the 
grey shaded region in Figure [3T|b). Such a method would be equivalent to the trapezoidal rule given in p.l| l. 
However, in our setting we would have to ignore, or subtract off, the area akeady accounted for in Region 
3, which is depicted as the shaded green section of Figure |3T|b). In doing so, the random variable needed 
in order to perform this step would necessarily be dependent upon the past (via Region 3), and our current 
analysis would break down. However, noting that Region 3 has the same area as Region 4, as depicted by the 
blue shaded region in Figure [TTJc), we see that it would be reasonable to expect that if one only uses Region 
5, as depicted as the grey shaded region in Figure [TTJc), then the accuracy of the method should be improved 
as we have performed a trapezoidal type approximation. Because 

(Region 5) ^ iv(o, (^^(^(O)) + 2^) |) ^ ^aUX{Q)) + 2V^N{Q,1) 

= ^2al{y*)-al{X{0))^N{0, 1), 

where V = (^\{y*) — "■^(-'^(0)), we see that this is precisely what is carried out by Step 2 of the Weak 
6'-Midpoint Trapezoidal Algorithm. 

3.2. A second point of view To obtain a higher order method one must both approximate well the 
expected drift term as well as the quadratic variation of the process. The basic idea of the Weak Trapezoidal 
Algorithm is to make a preliminary step using an Euler approximation and then use this step to make a higher 
order approximation to the drift integral and to the quadratic variation integral. Similar to ( |3.1| i the desired one 
step approximation to the quadratic variation integrals are 

/ al{X{s))ds^h 
Jo 



29 



2/ * 



5 



where all notation is as before. 

Considering just the variance terms of the quadratic variation, we let {e^} be an orthonormal basis and see 
that our method yields the approximation 



M 



Var(X(/;,) • eO « ^ Var(afc(yo)(i^fc • eO?7i/c + V ["i^K^*) " «2CTi(lo)] ^(^^fc • e,)jnkV ~ 0)h 
fc=i 

M 

fc=i 

If the step-size is sufficiently small then, [aiCT^ {y* ) — Q;20'^ (^o)] ^ is positive with high probability because of 
our uniform ellipticity assumption; and hence, 

M 
k=l 

which is a locally third order approximation to the true quadratic variation integral of 

lo 



Wai{X{h) ■ e,) = E V(i.fe • / al{X{s))ds. 



Notice that it was important in this simple analysis that the direction of variation ly^ stayed constant over the 
interval so that the two terms could combine exactly. Of course, one should really be computing the full 
quadratic variation, including terms such as Co\{X(h) ■ ei,X{h) ■ ej), but they follow the same pattern as 
above because each is a linear combination of the integral terms jj* af, {X{s))ds. 

4. Proof of Local Error Estimate We now give the proof of the local error estimate given in Theorem 



2.3 which is the central result of this paper. 

Proof, (of Theorem 2.3 1 We need to show that there exists a constant K so that for any f ^ one has 

\Ef{Yi) -Ef{X{h))\ < K\\f\\eh\ 

Hence for the reminder of the proof we fix an arbitrary f E C^. Observe that Step 1 of the Weak Trapezoidal 
Algorithm produces a value, y*, that is distributionally equivalent to y{9h), where y{t) solves 

M 

dy{t)^b{y{0))dt + J2<^k{ymiykdWt{t), 2/(0) = xq. (4.1) 
fc=i 

Likewise, Step 2 of the Weak Trapezoidal Algorithm produces a value, Yi, that is distributionally equivalent to 
y{h), where y{t) solves 

M 

dy{t) = («i6(2;*)-a26(xo))dt + 5]y^[aia2.(r)-a2^^(a:^o)]+ VkdWk{t), y{9h) = y* . (4.2) 

k=l 

Let J^t denote the filtration generated by the Weiner processes Wk{t) in ( |4. 1[ ) and ( |4.2| l. Then, 

Efiyih)) = E [E[f{yih)) \ Teh] ] =E [Eehf{y{h))], (4.3) 

where we have made the definition Eehl • ] =:E[ • | J-eh\- 
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Let A denote the generator for the process Bi denote the generator for the process ( |4.1| i, and B2 

denote the generator for the process (|4.2|i conditioned upon J^g/, . Then 



{Af)ix) = f[b]{x) + ^5]a^./"[z.,,i.,](x) 
fc 

(Bi/) (x) = f [6(xo)] {^) + \Y. '^fc (^0) V'l^^fe , i^^c] {x) 

k 

{B2f){x) = /'[ai6(2/*) - a2b{x^)]{x) + \ 5][ai(7fe(y*)2 - a^a^^x^fY r'\v^,M{x). 

k 

where f'[£,]{z) is the derivative of / in the direction evaluated at the point z. Note that {Af){xo) — 
{Bif){xo). For any integer fc > 2 we define recur sively {A''f){x) = {A{A''-^ f)){x), and similai'ly for 
Bi and B2- By repeated application of the Ito-Dynkin formula, see ifTTl . we have 

^ehf{y{h)) = f{y*) + [ E9h{B2f){y{s)) ds 
Jeh 

= f{y*) + {B2f){y*){l'0)h+ r (\eh{Blf){y{r)) dr ds 

J Oh J Oh 

= /(y*) + (S2/)(?/)(i - 9)h + {BlfW)^^^^ 

/ / ¥.eh{Blf){y{u))dudrds. 
eh Jeh Jeh 

The term {B2f)(y(u)) depends on the first six derivatives of /. Therefore, since f E 

i-h 



(4.4) 



V.eh{B'2f){y{u))dudrds 
eh Jeh Jeh 



< C\\f\\eh\ (4.5) 



for some constant C. Combining ( |4.3| l, ( |4.4| l, ( |4.5| l, and recalling that E/(Yi) = Ef{y{h)) gives 

Ef{Y,)=E[Eghfiyih))] = E/(y*) +E (i?2/)(y*)(l - e)h + E iBlf)iy*)^-^^^^ + 0{h'). (4.6) 

Here and in the sequel, we will write F = G + 0{hP) to mean that there exist a constant K depending on only 
a and b so that for all initial conditions xo 

\F-G\<K\\f\\ehP, (4.7) 

for h sufficiently small. In the spirit of the preceding calculation, repeated application of the Ito-Dynkin formula 
to produces 

EfiXih)) = fixo) + {Af){x,)h + (A2/)K)y + Oih"). 



The proof of the theorem is then completed by Lemma [4T| given below. Its proof, which is straightforward but 
tedious, is given in the appendix. □ 

Lemma 4. L Under the assumptions of Theorem 2.3 for all h > sufficiently small and f € one has 

E 



f{y*) + {B2fW){l - e)h + (L_^ 



--f{xo) + {Af){x,) + {A'f){xo)^ + 0{h^) . 



Remark 4.2. Comparing equation \'i.2\ and Lemma 4.1 shows that our algorithm can be viewed as providing 
an approximation to the two step Taylor series approximation. 
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5. Examples We present two examples that demonstrate the rate of convergence of the Weak Trapezoidal 
Algorithm with 6 — 1/2. In each example we shall compare the accuracy of the proposed algorithm to that of 
Euler's method and a "midpoint drift" algorithm defined via repetition of the following steps 



M 



Y, = + b{y*)h + ^ (Tfe(r,_i)z/fe % V^, 



(5.1) 



k=l 



where the notation is as before. We compare the proposed algorithm to that given via \5.\\ to point out that the 
gain in efficiency being demonstrated is not solely due to the fact that we are getting better approximations to 
the drift terms, but also because of the superior approximation of the diffusion terms. 



5.1. First Example. Consider the system 



dXi{t) 
dX2{t) 







dWi{t) 



1 

To 



dW2{t), 



(5.2) 



where Wi{t) and W2{t) are standard Weiner processes. In our notation bi{x) = xi, 62 (x) = 0, (Ti{x) = xi, 
o'2{x) — 1/10, and i^i = [0, l]'^, V2 = [1, 1]^. Note that the noise does not commute. It is an exercise to show 
that 



EX2(t)2 =EX2(0)2 - iEXi(0)2 + ^e2*(200EXi(0)^ 



1) 



t 

200 



1 

400' 



(5.3) 



For both Euler's method and the midpoint drift method \5.\\ we used step sizes hk = 1/3''', k e {1, 2, 3, 4, 5} 
and initial condition Xi{0) — ^2(0) = 1 to generate 500,000 sample paths of the system ( |5.2| i. We then 
computed 



errorfe(i) = EX2(t)2 



1 



5x10=" 



5 X 105 



(5.4) 



1=1 



where X '° (t) is the sample path generated numerically and EX2(i)^ is given via ( |5.3| l. We also generated 
10,000,000 sample paths using the Weak Trapezoidal Algorithm with the same initial condition and step 
sizes hk ~ l/(4fc), k E {1,2,3,4}. We then computed eiTorfe(t) similarly to before. The outcome of the 
numerical experiment is summarized in Figure 5.1a where we have plotted log(/ife) versus log(|errorfc(l)|) for 



the different algorithms. As expected, we see that the Weak Trapezoidal Algorithm gives an error that decreases 
proportional to h^, whereas the other two algorithms give errors that decrease proportional to h. 



5.2. Second Example. Now consider the following system that is similar to one considered in 11201 

dWiit) 
dW2{t), 



dXi{t) ' 




■ -X2{t) ' 


/sin2(Xi(i)H 


hX2(0)+6 


dX2{t) 




X,{t) 


+ V H 


- 1 








^ /cOs2(Xi(t)- 


f ^2(<))+6 








+ V t 


hi 



(5.5) 



where Wi{t) and W2{t) are independent Weiner processes. It is simple to show that 



E|X(t)p =EX(0)2 + 131og(l + t). 
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(5.6) 




-6 -5 -4 -3 -2 -1 -3 -2.5 -2 -1.5 -1 -0.5 



log h log h 

(a) First Example (b) Second Example 

Fig. 5.1: Log-log plots of the step-size versus the error for the three different algorithms. In (a) the example 
i5.2\ is considered. The best fit lines for the data (shown) have slopes 2.029, .998, and 1.030, for the Weak 
Trapezoidal Algorithm, Euler's method, and the midpoint drift method, respectively. In (b) the example in 
( |5.5| l is considered. The best fit lines for the data (shown) have slopes 2.223, .952, and 1.098, for the Weak 
Trapezoidal Algorithm, Euler's method, and the midpoint drift method, respectively. In both examples all 
results agree with what was expected. 



We used step sizes hk = l/(2fc), k G {1, 2, . . . , 8}, to generate five million approximate sample paths of the 
system ( |5.5| l using each of: (a) Weak Trapezoidal Algorithm, (b) Euler's method, and (c) the midpoint drift 
method ( |5.1[ l. We then computed 

error,(0=E|X(i)p-^^^ J2 

i—l 

where X^^ [t) is the sample path generated numerically and E|J'('(i)p is given via (|5.6[). The outcome is 



summarized in Figure 5.1b where we have plotted \og{hk) versus log(|errorfc(l)|) for the different algorithms. 
As before, we see that the Weak Trapezoidal Algorithm gives an error that decreases proportional to K^, whereas 
the other two algorithms give errors that decrease proportional to h. 

Remark 5.1. We note that in both examples we needed to average over an extremely large number of computed 
sample paths in order to estimate error^lt) for the Weak Trapezoidal Algorithm. This is due to the fact that the 
increased accuracy of the method quickly makes sampling error the dominant error 

6. The effect of varying 6 The term [aial{y*) — a20'^(li_i)] ^ in Step 2 of the Weak Trapezoidal 
Algorithm will yield zero, and the given step will have a local error of only O(h^), if 

a.aliy*) < a^aliY^-i) ^ aliv*) < ^al{Y,^^) = (1 - 20 + 20>2(r,_i)- 

We will call such a step a "degenerate" step. The function g{9) = 1 — 26* + 26^ is minimized at g{l/2) = 1/2, 



and g{6) ^ 1 as — or 6* 1. Therefore, as mentioned Remark 2.2 one would expect that as 6* ^ or 
6^1 more steps will be degenerate, and a decrease in accuracy, together with a bias against crfc decreasing, 
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0.2 0.4 0.6 0.8 1 ~ -3 -2.5 -2 -1.5 -1 

e log h 



(a) % of degenerate steps vs (b) accuracy for different 9 

Fig. 6.1: (a) The number of degenerate steps for the Weak Trapezoidal Algorithm applied to \6.\\ with h = 
1/10 and different values of 9. (b) The log h vs log(|error|) plot is given for different choices of 6 for the Weak 
Trapezoidal Algorithm applied to ( |5.2| i where the error is defined similarly to ( |5.4[ ). The best fit lines for the 
data (shown) have slopes 1.865, 1.996, 2.029, and 2.033 for 9 .05, .25, .50, .75, respectively. 



would follow. Using a step-size of h = 1/10, we tracked the percentage of degenerate steps for the simple 
system 



dX{t) = y/X{tY + 1 dW{t), X(0) = 1, (6.1) 
where W{t) is a standard Weiner process. To do so, we computed 10, 000 sample paths over the time interval 



[0, 1] for each of 6* = .02fc, k £ {1, . . . , 49}. The results are shown in Figure 6.1a where the behavior 
predicted above is seen. Curiously, the minimum number of rejections takes place at 9 = .42. It is also worth 
noting that one can check on computer software that in the general case the coefficient of for the one-step 
error grows like 1/6* as 6* ^ 0. This does not happen in the deterministic case ( |3.1| l. 

While the above considerations give some interesting insight into the effect of various 9, the situation is 
more complex. A 9 closer to one should give the method more stability, albeit at an expense as the rejection 
fraction increases as 9 approaches one. It would be interesting to perform a stability analysis in the spirit of IHl 



to better understand the effect of 9. In lieu of this. Figure 6. lb gives the result of a convergence analysis of the 
Weak Trapezoidal Algorithm applied to ( |5.2| l with different choices of 9. Interestingly, larger 9 seem to result 
in smaller (and hence better) convergence rate prefactors. This seems to indicate that in at least this example 
stability is an issue. 

The performance of the Weak Trapezoidal Algorithm as a function of ^^ is a topic deserving further consid- 
eration, but combining the above shows that 9 = 1/2 is a reasonable first choice, though stability considerations 
might lead one to consider a 9 closer to 1 . 

7. Comparison to Richardson Extrapolation It is illustrative to compare the Weak Trapezoidal Algo- 
rithm to Richardson extrapolation, which from a certain point of view is the method in the literature that is 
most similar to ours. See 1201 for complete details of Richardson extrapolation in the SDE setting. 

Let Zii/2{t) and Zh{t) denote approximate sample paths of ( |l.l| i generated using Euler's method with 
step sizes of h/2 and h, respectively. For all / satisfying mild assumptions, both Ef{Zi^/2{t)) and E/(Z/j(i)) 
will approximate E/(X(t)) with an order of 0{h). However, Richardson extrapolation may be used and the 
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hnear combination 2E/(Z,j /2(t)) - E/(Z/i) will approximate Ef{X{t)) with an order of 0(/i^) (see |!20|! ). 
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Fig. 7.1: The areas of space-time utilized by 2Z/j/2 ~ ^/i ™d the Weak Trapezoidal Algorithm for a single k 
and a single step. In 7.1 a), (T^(X(i)) increases and IZy^j^ — uses 77^^ + + 277^3, whereas the Weak 
Trapezoidal Algorithm uses r)Ai + + + '?A4- In the case when a\{X(i)) decreases, 7.1 ^b) above, the 
processes use 77^^ + r\A^ + JMs ^ '?A4 and r\Ax + ^7^2^ respectively. In both cases, it is the better use of the 
areas by the Weak Trapezoidal Algorithm that achieves a higher order of convergence. 



Of course, taking / to be the identity shows that the linear combination 2Z/i/2(i) — ^h(i) gives an 0{h^) 
approximate of the mean of the process. As Richardson extrapolation does not compute a single path, but 
instead uses the statistics from two to achieve a higher order of approximation for a given statistic, we will 
compare one step of the Weak Trapezoidal Algorithm with a step-size of h, to one step of size h of the process 
2Zi^/2{t) — Zh{t) with the clear understanding that 2Z;i/2(i)~-^?i(i) is only 0(/i) accurate for higher moments. 
Recall that systems of the form jl.l) are equivalent to those driven by space-time white noise processes 



p.3| l. As in Section 3.1 we consider how each method uses the areas of [0,cx))^ associated to Yk{du x ds) 
from p.3[ ) during one step. We will proceed considering a single k since it is sufficient to illustrate the point. 
For Ai C [0, c»)^, we denote by r/^. a normal random variable with mean and variance aiea{Ai). Recall that 
ijAi and rjAj are independent as long as Ai n Aj has Lebesgue measure zero. Consider ( |7.1| i(a) in which we are 
supposing that al{X{t)) increases over a single time-step. The change in the process due to this k would 
be i^k times 

VAi +VA2 +VA3- 

Similarly, the change in Zh would be i^k times rjAi +77A2 ■ Therefore, the change in the process 2Zi^/2 (t) — Z^ {t) 
would be Vk times 

VAi + VA2 + "^VAi ■ 

On the other hand, the change in the process generated by the Weak Trapezoidal Algorithm due to this k is 
times 

VAi +VA2 +VA3 +VAi- 

Therefore, and as expected, the means should be the same, but the variances should not as 



Var{2f]A3) = 4:Var{riA:,) = 2Var{TjA, +VaJ- 
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Similarly, in the case in which cr^(X(t)) decreases as depicted in ( |7.1| l(b), the process 2Zii/2{t) — Zh{t) would 
use r]ji^ + TjA^ + 77^3 — riA_^, whereas the Weak Trapezoidal Algorithm would use 77^^ + 77^^ . Again, the means 
will be the same, but the variances will not. In both cases, the Weak Trapezoidal Algorithm makes better use 
of the areas to approximate the quadratic variation of the true process, and thus achieves a higher order of 
convergence. 

8. Extension to General Uniformly Elliptic Systems For a moment let us consider the setting of general 
uniformly elliptic SDEs 

M 

dX{t) - h{X{t))dt + Y.9k{X{t)) dWu{t), 

k=i ^^-^l 

X{Q) = a; e M"* 

where 6 and are as before and : M'' M'' is such that if G(a:;) = {gi{x),--- , gM{x)){gi{x), ■ ■ ■ ,.gM(a;))^| 
then there exist positive A_ and such that 

< G{x)^ ■ ^ < x+\e 

for all 2;, ^ €E M''. For such a family of uniformly elliptic matrices a lemma of Motzkin and Wasow fT5\, whose 
precise formulation we take form Kurtz 1 10|, states that if the entries of G are C'^ then there exists an M and 
{CTfc : M'^ M>o : fc = 1, . . . , A/}, {i/fc e E'^ : fc = 1, • • • , Af} with (Jk e and strictly positive so that 

Hence ( |8.1| l has the same law on path space as with these ak and Vk- Of course M might be arbitrarily 
large (depending on the ratio of A+/A_) and hence it is more subtle to compare the total work for our method 
with a standard scheme based directly on ( |8.1[ ). Furthermore, depending on the dependence on x, it is not 
transparent how to obtain the vectors v and functions a exactly. Approximations could be obtained using the 
SVN of the matrix G{x) for fixed x but we do not explore this further here. 

9. Conclusions and Further Extensions We have presented a relatively simple method directly appli- 
cable to a wide class of systems which is weakly second order. We have also shown how, at least theoretically, 
it should be applicable to systems which do not satisfy our structural assumptions but are uniformly elliptic. 
We have picked a particularly simple setting to perform our analysis to make the central points clearer. The 
assumption that h and ak are uniformly bounded can be relaxed to a local Lipschitz condition. That is to say, if 
b and a and their needed derivatives are not bounded uniformly, but rather are bounded by an appropriate Lya- 
punov function, then it should be possible to extend the method directly to the setting of unbounded coefficients 
provided the method is stable for the given SDE (see for instance [12|). If the SDH is not globally Lipschitz 
then using an implicit drift split-step method as in [12J, an adaptive method as in 1 11 1, or a truncation method 
as in lfT4l should extend to our current setting. More interesting is relaxing the non-degeneracy assumption on 
the CTfe, which was used to minimize the probability of the diffusion correction being negative. This tact is in 
some ways reminiscent of 1 14| in that a modification of the method is made on a small set of paths, though the 
take here is quite different. It would be interesting to use the probability that the correction to the diffusion is 
negative to adapt the step-size much in the spirit of |11 1. Lastly, there is some similarity of our method with 
predictor corrector methods. In the deterministic setting, predictor corrector methods not only have a higher 
order of accuracy but also have better stability properties. There have been a number of papers exploring this 
in the stochastic context (see IS] IH [191 El). It would be interesting to do the same with the method presented 
here. 
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Appendix A. Proof of Lemma [43| The proof of Lemma [4T| requires the replacement of the terms of the 

form [Q;icr|(y*) — a2cr|(a;o)]^ with [aial{y*) — a20'f(xo)]. The following two lemmas show that this can be 
done at the cost of an error whose size is 0{h^). Here 0{h^) has the same meaning described earlier around 
We begin with an abstract technical lemma where p and q satisfy 1/p + 1/q = 1. 



Proposition A.l. Let X and Y be a real valued random variables on a probability space with 
\XY\LPl^n) < oo for some p e (l,oo]. Then \EY[X]+ -EYX\ < \YX\Lp(^n){V{X < 0})i/9. Similarly if 
X,Y and Z are real valued random variables with \ZXY\i^p(ji) < oo and A = {X < 0} U {Z < 0} then 

\EY[X]+[Z]+ -EYXZ\ < 2|ZyX|ip(o)(P{v4})i/9. 

Proof Let A = {X < 0} and 9 = p/{p - 1). Then |Ey([X]+ - X)\ < E\Y\\[X]+ - X\1a < 
|rX|ip(o)(P(^))^/«, showing the first claim. For the second notice that Ey[X] + [Z] + -ErXZ = (Er[X] + Z-| 
EYXZ) + (EF[X] + - EF + and that each of the terms in parentheses can be bounded by the first 
result. □ 

Corollary A. 2. Let G with infx crk{x) > Q for all k and let Y be a random variable with \Y\ < C 
a.s. for some C. Then for any p > 1 there exists an Hq so that 

EY[a,al{y*)-a2alixo)]+ - EY[ai<jl(y*) - a2<jl{xo)] + 0{hP) 

EY[aialiy*)-a2al{xo)] + [aiC7j{y*) ~ a2CT|(a;o)]+ = EY[aial{y*) - a2al{xo)][aia]{y*) - azaf (xq)] + 0(/iP)| 
for all h G (0, Hq] and k,£ £ {!...., AI}, where y* is defined via Step 1 of the Weak Trapezoidal Algorithm. 

012 I I 

Proof. Define the event Ak — {ckiy*) < — crkixo)}- In light of Proposition A.l it is sufficient to show 

ai I I 

that for any p > 1 there exists a Cp such that V{Ak) < Cph^. Because ak is Lipschitz there exists a positive C 

such that 

alixo +S)- ^alixo) > (1 - "')a^(xo) - C|<5|, 
ai ai 

for any (5 > 0. In particular, setting S = y* — xq = b{xo)Oh + J2j f j (xq)\/ Oh Vj rii^j , and noting that a2 < ai 
and that the cr's are uniformly bounded from both above and below, the result follows from the Gaussian tails 
of the 77's. □ 



Proof, (of Lemma 4. 1 1 From Taylor's theorem and the definition of the operators involved one has 



Kfiy*) = f{xo) + iBJ){xo)eh + {Blf){x^f^ + 0{h^) 
= /(xo) + {Af){x,)6h + {Blf){x/^ + 0{h^) . 

In the last line, we have used the observation that (i?i/)(xo) — (A/)(a;o). Now we turn to E{B2f){y*). We 



begin by using Lemma A. 2 to remove the [ • ]+. Then we use the fact that ai — 0:2 = 1 and Taylor's theorem 
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to expand vaiious terms to produce the following: 

E{B2f){y*) = Ef{y*)[aMy*) - a^bix,)] + ^EY,M{y*) ~ a-,al{x^)Y f'X^k. ^kW) 

k 

= E/'(y*)[ai5(y*) - a^bix^)] + ^Ej2[a,<jl{y*) - a2al{xo)]f"Wk,>^k]{y*) + 0{h') 

k 

fc 

+ EBi - a^bixo)] + ^ ^^("I'^fc ~ "2'T^(a:o))/"K, ^^fc]) {xo)9h + 0{h^) 

k 

= {Af){xo) + ai{BMmxo)eh - a2{Blf){xo)eh + 0{h^) 
= {Af)ixo) + a,{A^f){xo)eh - a2{Bff)(xo)0h + 0{h^). 

Similar reasoning produces 

nBlf){y*) = E(i?2(/'[«ife(y*) - a^bixo)] + \ Y^o^i^liv*) " a2ol{x^)]+ f"[vu,Vk]){y*)) 

k 

- r[b{xo),b{xo)]{xo)+EY,[c^i<^l{y*) - a2<jl{xo)]+r'Wk,'^k,b{xo)]{xo) 

k 

+ ^E^[aicr^(y*) - a2CTfc(xo)] + [ai(T|(2/*) - a2(T'^j{xo)\^ f""[vu,Vk,Vj,Vj]{xo) + 0{h) 
= f'Hxo), b{xQ)]{xo) + ^ cr^(xo)/"'K, i^fc, b{xo)]{xo) 

k 

+ \ ^'^k{xo)<^'^(xo)f'"Wk,i^k,i'j,i^j]{xo) + 0{h) 
= iBlf){xo) + 0{h). 

Combining these estimate and the fact that 2(1 — 9)6a2 = 9"^ + {1 — O)"^ and 2(1 — 0)6ai = 1, produces the 
quoted result after some algebra. □ 

Appendix B. Operator Bound for Vt- ^ C'^. 

In this section, we show that if 6, ai £ then Vt is a bounded operator from C™ to C™ for m £ 
{0, • • • , fc}. The fc = case follows immediately from |/(a:)| < ||/||o for all x E M'^. To address the higher k, 
we introduce the first k variations of equation ( |l.l| i. 

For any ^ e M'' we denote the first variation of ( |l.l| i in the direction ^ by J'^^\t, x)[^] which solves the 
linear equation 

M 

dj^') it, x) [i] - (v&) {x(t) ) [ (<, X) m dt+Y. (va,) (x(i)) [ (<, x) m dWk {t) , 

fc=i 

/i)(0,a;)[C] = C and X{0) ^ x 

Similarly for ^ = (^1,^2) € the second variation of X{t) (in the directions ^ will be denoted by 
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J(2)(i,x)[^] and defined by 



M 



fc=i 

M 



{y%){x{t))[j^'\t,x)[^{\,j'^^\t,x)[^] + Y.^y^''k){xmj^^^ 



k=l 



j(^)(0,a;)[C] = and X{0) ^ x . | 

These equations were obtained from successive formal differentiation of By further formal differentiation 
we obtain analogous equations for the fc-variation x)[^] where — (^i, • • • ,£,k) € K'^ is the vector of 

directions. It is a standard fact that if the coefficients &, aj are in C*^ then for any i > 

supE,sup{ sup |j(")(s,a;)[a,...,Cn]r:C»eK'^with|$,| = l} <oo. 
X se[o,t] 



This can be found in Lemma 2 in Q on p. 196 or in a slightly different context in Proposition 1.3 in |[T6] 
With these definitions in hand, we have that for any / S that 

V(PJ)(a;)K]=E,/'(X(t))[j(i)(t,x)K]], 

V'(PJ)(x)K]=E,/'(X(i))[j(2)(t,x)K]]+E,/(2)(x(i))[j(i)(<,x)Ki],j(ini,2:)[6]]. 
Using the moment bounds we have that for q > 1 and an ever changing constant C, 

E sup \V{Vtf){xm\'^ <C||/||^i sup < < oo 

l?l=i 

E sup |v2(7'J)(x)Ki,6]r <C||/||^.((E sup |jW(i,x)[a]|'')^+E sup 

<C||/||^. <oo 

Continuing in this manner we see that for any positive integer m if /, b, ui G C™ then for any g > 1 one has 

E sup iv"(7'*/)(x)[a,--- ,em]r <cii/rc„ <oo 
i?.i=i 

for some C . Now observe that taking q — \ proves the desired claim on the operator norm of Vt from to 
since 

fc k 

\\Vt!\\k<cy^ sup |(v^7't/)(x)[ei,--- <cVii/iic. <cmc- 

j=0 l«'l = i 
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